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Using holography, we study the collision of planar shock waves in strongly coupled N = 4 super- 
symmetric Yang-Mills theory. This requires the numerical solution of a dual gravitational initial 
value problem in asymptotically anti-de Sitter spacetime. 



Introduction. — The recognition that the quark-gluon 
plasma (QGP) produced in relativistic heavy ion col- 
lisions is strongly coupled [T], combined with the ad- 
vent of gauge/gravity duality (or "holography") |3j, 
has prompted much work exploring both equilibrium and 
non-equilibrium properties of strongly coupled Af = 4 su- 
persymmetric Yang-Mills theory (SYM), which may be 
viewed as a theoretically tractable toy model for real 
QGP. Multiple authors have discussed collisions of in- 
finitely extended planar shock waves in SYM, which may 
be viewed as instructive caricatures of collisions of large, 
highly Lorentz-contracted nuclei. In the dual description 
of strongly coupled (and large N c ) SYM, this becomes a 
problem of colliding gravitational shock waves in asymp- 
totically anti-de Sitter (AdSs) spacetime. Previous work 
has examined qualitative properties and trapped surfaces 
[IHZ], possible early time behavior [8HTU]. and expected 
late time asymptotics [ITJ [12] • As no analytic solution is 
known for this gravitational problem, solving the gravita- 
tional initial value problem numerically is the only way to 
obtain quantitative results which properly connect early 
and late time behavior. In this letter, we report the re- 
sults of such a calculation, and examine the evolution of 
the post-collision stress-energy tensor. 

Unlike previous work considering singular shocks with 
vanishing thickness [9], or shocks driven by non- 
vanishing sources in the bulk [5J [5] , we choose to study 
planar gravitational "shocks" which are regular, non- 
singular, source-less solutions to Einstein's equations. 
Equivalently, we study the evolution of initial states in 
SYM with finite energy density concentrated on two pla- 
nar sheets of finite thickness (and Gaussian profile) , prop- 
agating toward each other at the speed of light. The 
dual description only involves gravity in asymptotically 
AdSs spacetime; the complementary 5D internal mani- 
fold plays no role and may be ignored. Consequently, our 
results apply to all strongly coupled 4D conformal gauge 
theories with a pure gravitational dual description. 

Gravitational description. — Diffeomorphism invari- 
ance plus translation invariance in two spatial directions 
allows one to write the 5D spacetime metric in the form 

ds 2 = -Adv 2 + Y? [e B dx\ + er 2B ' dz 2 ] + 2dv (dr + Fdz), 

(1) 



where A, B, E, and F are all functions of the bulk radial 
coordinate r, time v, and longitudinal coordinate z. We 
employ generalized infalling Eddington-Finkelstein coor- 
dinates. Lines along which all coordinates except r are 
constant are infalling radial null geodesies; the radial co- 
ordinate r is an affine parameter along these geodesies. 
At the boundary, located at r — oo, v coincides with 
time in the dual quantum field theory. The geometry in 
the bulk at v > is the causal future of v — on the 
boundary. The ansatz ([I]) is invariant under the residual 
diffeomorphism r — > r + £, with £ an arbitrary function 
of v and z. 

For a metric of the form , Einstein's equations (with 
cosmological constant A = — 6) imply 

= S" + |(B') 2 S, (2a) 

= E 2 [F" - 2{d 3 B)' - ZB'd 3 B] + 4E'd 3 E , 

- S [3E'_F' + 4(<f 3 E)' + 6B'd 3 T,] , (2b) 

= E 4 [A" + SB'd+B + 4] - 12E 2 S'd + I] 
+ e 2B {£ 2 [\{F') 2 -l(d 3 B) 2 -2d 2 B\ 

+ 2(d 3 S) 2 - 4E [2(d 3 B)d 3 E + d 2 £] } , (2c) 

= 6S 3 (d+S)' + 12S 2 (S'd+E - E 2 ) - e 2B |2(d 3 E) 2 
+ Z 2 [±(F') 2 + (d 3 Fy+2F'd 3 B-l(d 3 B) 2 ~2d 2 B] 
+ E [(F'-8d 3 B) d 3 E - 4<i 2 E] } . (2d) 

= 6E 4 (d+B)' + 9E 3 (E'd + S + B'd + T) 

+ e 2B {Y?[{F') 2 +2{d 3 F)' +F 1 d 3 B-{d 3 B) 2 -d 2 3 B] 
+ 4(d 3 E) 2 - E [{AF'+d 3 B) d 3 E + 2d 2 E] } , (2e) 

= GE^E - 3E 2 A'd + E + 3E 3 (d + B) 2 
- e 2S { (d 3 E + 2Hd 3 B){2d + F + d 3 A) 

+ ^[2d 3 {d + F)+d 2 A]} 1 (2f) 

= E [2d+(d 3 E) + 2d 3 {d + H) + 3F'd+E] 

+ E 2 [d+(F') + d 3 {A r ) + 4d 3 (d+B) - 2d+(d 3 B)} 
+ 3E (Ed 3 S + 2d 3 E) d+B - 4(d 3 E)rf + E , (2g) 

where, for any function h(v, z, r), h! = d r h and 

d + h = d v h+ hAd r h, d 3 h = d z h — F d r h . (3) 
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Note that h! is a directional derivative along infalling ra- 
dial null geodesies, d+h is a derivative along outgoing 
radial null geodesies, and d$h is a derivative in the lon- 
gitudinal direction orthogonal to both radial geodesies. 

Near the boundary, Einstein's equations may be solved 
with a power series in r. Solutions with flat Minkowski 
boundary geometry have the form 



A = 



1+| + ^ + £+0(r-«) 
h 



F = d z i + i\ + 0(r" 3 ) . 

B= b ^ + 0(rS), 
r 4 

£ = r + £ + 0(r- 7 ) , 



(4a) 

(4b) 

(4c) 
(4d) 



The coefficient £ is a gauge dependent parameter which 
encodes the residual diffcomorphism invariance of the 
metric. The coefficients a 4 , b 4 and f 2 are sensitive to 
the entire bulk geometry, but must satisfy 



d v a 4 



d z f 2 , d v f 2 = -d z (\( H + 2b i ). (5) 



These coefficients contain the information which, under 
the holographic mapping of gauge/gravity duality, de- 
termines the field theory stress-energy tensor [T5] , 
Defining £ = ^T°°,V±= ^ T ±A -, S=^T°', and 



V\\ = ^-T zz , one finds 



N- 



£ = 
S = -h 



4^4 



Vx = -|a 4 + b 4 , 
"Pi, = -\a 4 -2b 4 . 



(6a) 
(6b) 



Eqs. Q and ^ imply d^T^ = and = 0. 

Numerics overview. — Our equations ([2]) have a natu- 
ral nested linear structure which is extremely helpful in 
solving for the fields and their time derivatives on each 
v = const, null slice. Given B, Eq. (2a I may be inte- 



grated in r to find E. With B and £ known, Eq. (2b) 



may be integrated to find F. With B, £ and F known, 
Eq. Q may be integrated to find d+E. With B, E, F 



and d+E known, Eq. (2e| may be integrated to find d+B. 
Last, with B, E, F, d + E and d+B known, Eq. (2c) may 



be integrated to find A. At this point, one can compute 
the field velocity d v B = d + B — \AB' , evolve B forward 
in time to the next time step, and repeat the process. 

In this scheme, each nested equation is a linear ODE 
for the field being determined, and may be integrated in 
r at fixed v and z. The requisite radial boundary condi- 
tions follow from the asymptotic expansions Q. Con- 
sequently, the initial data required to solve Einstein's 
equations consist of the function B plus the expansion 
coefficients a 4 and f 2 — all specified at some constant v 
— and the gauge parameter £ specified at all times. Val- 
ues of o 4 and / 2 on future time slices, needed as boundary 
conditions for the radial equations, are determined by in- 
tegrating the continuity relations ^ forward in time. 



Eqs. <[2f|) and ( |2g| are only needed when deriving 
the series expansions Q and the continuity conditions 
([5]). In this scheme, they are effectively implemented as 
boundary conditions. Indeed, the Bianchi identities im- 
ply that Eqs. 



(pB and (2g 



are boundary constraints; if 
they hold at one value of r then the other Einstein equa- 
tions guarantee that they hold at all values of r. 

An important practical matter is fixing the computa- 
tional domain in r. If an event horizon exists, then one 
may excise the geometry inside the horizon, as this re- 
gion is causally disconnected from the outside geometry. 
Moreover, one must excise the geometry to avoid singu- 
larities behind the horizon [Tl] . To perform the excision, 
we identify the location of an apparent horizon (an outer- 
most marginally trapped surface) which, if it exists, must 
lie inside an event horizon [TS] . For the initial conditions 
discussed in the next section, the apparent horizon al- 
ways exists — even before the collision — and has the 
topology of a plane. Hence, one may fix the residual dif- 
feomorphism invariance by requiring the apparent hori- 
zon position to lie at a fixed radial position, r = 1. The 
defining conditions for the apparent horizon then imply 
that fields at r = 1 must satisfy 



= 3E 2 rf+E -d z (FE e 2B ) + §F 2 E'e 



2L1 



(7) 



which is implemented as a boundary condition to deter- 
mine £ and its evolution. Horizon excision is performed 
by restricting the computational domain to r € [l,oo]. 

Another issue is the presence of a singular point at 
r = 00 in the equations ([2]). To handle this, we discretize 
Einstein's equations using pseudospectral methods |16j . 
We represent the radial dependence of all functions as a 
series in Chebyshev polynomials and the z-dependence 
as a Fourier series, so the z-direction is periodically com- 
pactified. With these basis functions, the computational 
domain may extend all the way to r = 00, where bound- 
ary conditions can be directly imposed. 

Initial data. — We want our initial data to describe two 
well-separated planar shocks, with finite thickness and 
energy density, moving toward each other. In Fefferman- 
Graham coordinates, an analytic solution describing a 
single planar shock moving in the =Fz direction may be 
easily found and reads [TT] . 

ds 2 = r 2 {-dx+dx- + dx 2 A + \ [dr 2 + h(x±)dsc±] , (8) 

with x± = t ± z, and h an arbitrary function which we 
choose to be a Gaussian with width w and amplitude /j, 3 , 



h(x ± ) = ^ [2nw 2 r 1 ' 2 e-^i^ 



(9) 



The energy per unit area of the shock is fj, 3 (N 2 / 2ir 2 ) . If 
the shock profile h has compact support, then a super- 
position of right and left moving shocks solves Einstein's 
equations at early times when the incoming shocks have 
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FIG. 1: Energy density £/fi 4 as a function of time v and 
longitudinal coordinate z. 

disjoint support. Although this is not exactly true for our 
Gaussian profiles, the residual error in Einstein's equa- 
tions is negligible when the separation of the incoming 
shocks is more than a few times the shock width. 

To find the initial data relevant for our metric ansatz 
([!]) , we solve (numerically) for the diffeomorphism trans- 
forming the single shock metric ([8| from Fcfferman- 
Graham to Eddington-Finkelstein coordinates. In par- 
ticular, we compute the anisotropy function B± for each 
shock and sum the result, B = B + + B-. We choose the 
initial time Vq so the incoming shocks are well separated 
and the B± negligibly overlap above the apparent hori- 
zon. The functions 04 and ji may be found analytically, 

a 4 = -| [h(v +z)+h(v -z)} , f 2 = h(v +z)-h(v -z). 

(10) 

A complication with this initial data is that the metric 
functions A and F become very large deep in the bulk, 
degrading convergence of their spectral representations. 
To ameliorate the problem, we slightly modify the initial 
data, subtracting from 04 a small positive constant 5. 
This introduces a small background energy density in 
the dual quantum theory. Increasing S causes the regions 
with rapid variations in the metric to be pushed inside 
the apparent horizon, out of the computational domain. 

We chose a width w — 0.75/ /i for our shocks. The 
initial separation of the shocks is Az = 6.2/ /t. We chose 
S = 0.014 /1 4 , corresponding to a background energy den- 
sity 50 times smaller than the peak energy density of the 
shocks. We evolve the system for a total time equal to 
the inverse of the temperature associated with the back- 
ground energy density, Tbkgd = 0.11//. 

Results and discussion. — Figure [I] shows the energy 
density £ as a function of time v and longitudinal position 
z. On the left, one sees two incoming shocks propagating 
toward each other at the speed of light. After the colli- 
sion, centered on v = 0, energy is deposited throughout 
the region between the two receding energy density max- 
ima. The energy density after the collision does not re- 
semble the superposition of two unmodified shocks, sepa- 
rating at the speed of light, plus small corrections. In par- 




FIG. 2: Energy flux S / fx as a function of time v and longi- 
tudinal coordinate z. 
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FIG. 3: Longitudinal and transverse pressure as a function 
of time v, at z = and z — 3//i. Also shown for compari- 
son are the pressures predicted by the viscous hydrodynamic 
constitutive relations. 

ticular, the two receding maxima are moving outwards at 
less than the speed of light. To elaborate on this point, 
Figure [2] shows a contour plot of the energy flux S for 
positive v and z. The dashed curve shows the location 
of the maximum of the energy flux. The inverse slope 
of this curve, equal to the outward speed of the maxi- 
mum, is V — 0.86 at late times. The solid line shows the 
point beyond which S/fi 4 < 10~ 4 , and has slope 1. Ev- 
idently, the leading disturbance from the collision moves 
outwards at the speed of light, but the maxima in £ and 
S move significantly slower. 

Figure[3]plots the transverse and longitudinal pressures 
at z — and z = 3//i, as a function of time. At z = 0, 
the pressures increase dramatically during the collision, 
resulting in a system which is very anisotropic and far 
from equilibrium. At v — — 0.23//i, where V\\ has its 
maximum, it is roughly 5 times larger than V±. At late 
times, the pressures asymptotically approach each other. 
At z = 3//i, the outgoing maximum in the energy density 
is located near v = There, V\\ is more than 3 times 
larger than 'P±. 

The fluid/gravity correspondence [T7] implies that at 
sufficiently late times the evolution of T^ v will be de- 
scribed by hydrodynamics. To test the validly of hydro- 
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dynamics, Fig.[3]also plots (as dashed lines) the pressures energy density. 



phydio an( j p j o p rec jj c ^ ec [ by ^ e first-order viscous 

hydrodynamic constitutive relations [TS]. At z = the 
hydrodynamic constitutive relations hold within 15% at 
time i>h y dro = 2.4//i, with improving accuracy thereafter. 
At z — 3/fi, they hold within 15% or better accuracy 
after v hy dro = 2.1/fi. 

At z = 0, where the flux 5 = 0, the constitutive rela- 
tions imply that the difference between 'P 1 y ro and r p^ ydro 
is purely due to viscous effects. Fig.[3]shows that there is 
a large difference between V± and Vn when hydrodynam- 
ics first becomes applicable, implying that viscous effects 
are substantial. 

We have also examined the influence of second-order 
corrections in the hydrodynamic constitutive relations. 
At z = 0, the second-order corrections only change Vhydro 
by about 1%, whereas at z = 3//i their addition increases 
^hydro by 20%. Evidently, in front of the receding maxima 
in £ and 5, second-order corrections are large, the system 
is still far from equilibrium, and agreement with the first- 
order constitutive relation is largely fortuitous. 

Hydrodynamic simulations of heavy ion collisions at 
RHIC suggest that the produced QGP thermalizes in a 
time perhaps shorter than lfm/c [19[. Crudely modeling 
the boosted gold nuclei by our translationally invariant 
Gaussian shocks, for RHIC energies we estimate /i ~ 2.3 
GeV. Our results from Fig. [3] then imply that the total 
time required for apparent thermalization, from when the 
Gaussian shocks start to overlap significantly to the onset 
of validity of hydrodynamics, is Ai) to t ~ 4//i ~ 0.35 fm/c. 
Similar results for far-from-equilibrium relaxation times 
in SYM were also reported in Ref . [2U] . 

We conclude by discussing the effect of the background 
energy density on our results. In the absence of a colli- 
sion, the presence of the background energy density will 
cause a propagating shock to slowly attenuate in ampli- 
tude and eventually thermalize. We have computed sin- 
gle shock propagation with the background energy den- 
sity used above, and found that the shock attenuates in 
amplitude by 2.5% in a time Av = 1/lbkgd- We have also 
studied the effect of increasing or decreasing the back- 
ground energy density by a factor of 1.5. This results in 
a spacetime dependent 0(S) change in the stress tensor 
and perturbs the thermalization time fhydro 

(at 2 = 0) 

by 1%. Lastly, we have also studied an expanding sheet 
of plasma which is initially localized at z = and sur- 
rounded by vacuum. The initial conditions consisted of 
B = 0, a Gaussian profile in the energy density, and 
vanishing energy flux. At late times the stress has two 
outward moving maxima, just as it does for the colliding 
shocks. Furthermore, the tails in the stress move out- 
wards at the speed of light whereas the maxima move at 
a speed around 15% less. Consequently, we are confident 
that the deviation of the speed of the maxima from 1 in 
Figs. [I] and [2] is not an artifact caused by the background 
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